/**
 * @file VC_1d.h
 * @brief The VC for one dimension
 * @author XDDDD
 * @version 
 * @date 2021-06-10
 */

#ifndef __PAISLEYPARK__VC_1D_H__
#define __PAISLEYPARK__VC_1D_H__


#include "MultigridSolver.h"

class VC_1d : public MultigridSolver<1>{
public:
	VC_1d() : MultigridSolver<1>(){};
	~VC_1d(){};
	VC_1d(std::string _RestrictionId, std::string _InerpolationId) : MultigridSolver<1>(_RestrictionId, _InerpolationId){};
	void Solve(VectorFunction<1>& _f, VectorFunction<1>& _g, std::vector<double>& _v,
			const size_t _N, const size_t _M, const size_t itr_num, 
			const size_t nu1, const size_t nu2, const double _epsilon) override;
};


void VC_1d::Solve(VectorFunction<1>& _f, VectorFunction<1>& _g, std::vector<double>& _v, const size_t _N,
				  const size_t _M, const size_t itr_num, const size_t _nu1, const size_t _nu2, const double _epsilon) {
	bool is_stop = false;
	this->set(_v, _N, _f);
	size_t N = this->__U.get_N();
	size_t pre_N = N;
	for(size_t j = 0; j != itr_num; j++) {
		is_stop = false;
		while(!is_stop) {
			for(size_t i = 0; i != _nu1; i++) {
					this->smooth(pre_N);
			}
			if(pre_N == _M) {
				is_stop = true;
				break;
			}
			this->set_e(pre_N);
			
			pre_N -= 1;
			this->__U.clean(pre_N);
		}
		this->solve_coarest(pre_N);
		is_stop = false;
		while(!is_stop) {
			this->putback_e(pre_N);
			pre_N++;
			for(size_t i = 0; i != _nu2; i++) {
				this->smooth(pre_N);
			}
			if(pre_N == N) {
				is_stop = true;
			}
		}
	}

	std::vector<double> x(1, 0);
	std::vector<size_t> coord(1, 1);
	size_t pole = 0;
	for(size_t i = 1; i != (1<<_N); i++) {
		x[0] = 0;
		__U(i, _N) += (1.0 - 1.0*i/(1 << _N))*_g(x);
		x[0] = 1;
		__U(i, _N) += 1.0*i/(1 << _N)*_g(x);
	}

};

#else
//Do nothing.
#endif
